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Abstract 

Non-ideal fluid dynamics with cylindrical symmetry in transverse direction and 
longitudinal scaling flow is employed to simulate the space-time evolution of the 
quark-gluon plasma produced in heavy-ion collisions at RHIC energies. The dy- 
namical expansion is studied as a function of initial energy density and initial time. 
A causal theory of dissipative fluid dynamics is used instead of the standard theo- 
ries which are acausal. We compute the parton momentum spectra and HBT radii 
from two-particle correlation functions. We find that, in non-ideal fluid dynamics, 
the reduction of the longitudinal pressure due to viscous effects leads to an increase 
of transverse flow and a decrease of the ratio R ou t/R s ide a s compared to the ideal 
fluid approximation. 
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1 Introduction 



High-energy heavy-ion collisions offer the opportunity to study hot and dense nuclear 
matter [1]. Lattice QCD calculations [2] of thermodynamical functions indicate that, 
at zero net baryon number density, strongly interacting matter undergoes a rapid 
transition from (a chiraliy broken, confined) hadronic phase to a (chirally symmetric, 
deconfined) quark-gluon plasma (QGP). The transition temperature is found to be 
around T c ~ 160 MeV. One of the primary goals of relativistic heavy-ion physics is 
the creation and observation of the predicted QGP. 

Fluid dynamics is the only dynamical model which provides a direct link between 
collective observables and the equation of state. Besides, it can also be linked to mi- 
croscopic models via the energy-momentum tensor or via the transport coefficients. 
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The applicability of ideal fluid dynamics rests on the assumption of local thermo- 
dynamic equilibrium. That is, in the local rest frame of the fluid the single-particle 
phase space distribution function is assumed to be isotropic in momentum space. 
However, for the early stage of a collision this may not be a very accurate ap- 
proximation because the initial particle momenta are predominantly oriented in the 
longitudinal (beam) direction. Then it is necessary to account for the anisotropy of 
the phase space distribution in order to carry out a detailed analysis of the dynamics 
of the matter. Because of the anisotropy, shear viscous stresses arise. In the case of 
strongly interacting mixtures, in addition to shear and bulk viscous pressure, also 
diffusion has to be included in the analysis [3,4]. 

While for ideal fluid dynamics we only consider the conservation equations, in non- 
ideal fluid dynamics we need to include a second set of transport equations. The 
standard theories of irreversible processes due to Eckart [5] or Landau [6] exhibit 
undesirable effects (i.e., they predict an infinite speed of propagation for thermal 
and viscous signals and unstable equilibrium states [7]). It is then necessary to re- 
sort to another non-ideal fluid-dynamical theory that does not have this anomalous 
behavior. Causal theories for non-ideal fluids were formulated by Miiller [8] based on 
Grad's 14-moment method [9] and written in relativistically covariant form by Is- 
rael and Stewart [10]. These are based upon extended irreversible thermodynamics, 
where the entropy current vector of standard thermodynamics is extended by in- 
cluding terms that are quadratic in the dissipative quantities. Causal fluid dynamics 
seems to be a good candidate to replace the standard Eckart or Landau theory. It 
has already been applied previously in Ref. [11] to one-dimensional Bjorken scaling 
flow [12]. 

In the present work we extend the investigation of Ref. [11] to more realistic, three- 
dimensional geometries. We consider a system consisting of quarks and gluons with 
an equation of state of an ultrarelativistic ideal gas. This implies that the bulk vis- 
cosity vanishes [13]. We also assume the system to be net baryon free. Consequently 
the heat flow is small and will be neglected. The dominant dissipative corrections are 
then due to shear viscosity. The evolution of the system is studied assuming trans- 
verse cylindrical symmetry and Bjorken scaling flow in the longitudinal direction. 
To explore observable consequences of non-ideal fluid behaviour, we calculate the 
transverse momentum distribution of the produced partons and the two-particle cor- 
relation functions in side-, out- and long-direction, giving rise to the so-called HBT 
radii R S ide, Rout, and Ri ong . At this point, we do not consider the hadronization 
of partons since the transport coefficients in the mixed phase are unknown. Nev- 
ertheless, the present work is the first to solve the second order (and thus causal) 
equations of relativistic dissipative fluid dynamics for a realistic three-dimensional 
geometry. Previous attempts use the first order equations which are acausal as a 
matter of principle. 

The remainder of the paper is organized as follows. In section 2 we briefly review 
the equations of causal non-ideal fluid dynamics in a general frame. In section 3 
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we present the relaxational transport equations necessary to solve non-ideal fluid 
dynamics (also in a general frame). In section 4 we explicitly present the equations 
for the Bjorken cylinder-type expansion. In section 5 we outline the initial conditions, 
equations of state and transport coefficients used in the calculations. In section 
6 we discuss the freeze-out prescription for particle spectra and HBT radii. We 
also discuss the dissipative corrections to the distribution function. In section 7 we 
present our results for the spectra and HBT radii. In section 8 we summarize our 
work. Mathematical details are deferred to the Appendix. 

Our units are h — c — = 1- The metric tensor is g^ v = diag(+, — , — , — ). The 
scalar product of two 4- vectors a M , 6 M is denoted by a M g^ v b u = 6 M , and the scalar 
product of two 3- vectors a and b by a ■ b. The notations = (A a/3 + A /3a )/2 and 
A\ a fl = (A a/3 — A l3a )/2 denote symmetrization and antisymmetrization, respectively. 
The 4-derivative is denoted by d a = d/dx a . An overdot denotes A = u x d\A. 



2 Basics of non-ideal fluid dynamics 

In this section we present the equations of causal non-ideal fluid dynamics. We 
consider a fluid that consists of a single conserved charge. For a mixture of several 
conserved charges see Ref. [4]. The variables of concern are the (net) charge current 
iV M , the energy-momentum tensor T^ v and the entropy flux S^. The divergence of 
the energy- momentum tensor and the charge current vanish, i.e., energy-momentum 
and charge are conserved. However, in general the divergence of the entropy current 
does not vanish. The second law of thermodynamics requires that it be a positive 
and nondecreasing function of time. In equilibrium the entropy is maximum and the 
divergence of the entropy current vanishes. Thus, 



d v T' u ' = Q , 
d^>0 , 



(1) 
(2) 
(3) 



where 




(4) 
(5) 



W v W vA 



w 



) 
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In the local rest frame defined by w M = (1,0) the quantities appearing in the above 
tensors have the following meaning. Defining the projection operator A^ u = g^ u — 
u^u v = A Ufl onto the 3-space orthogonal to w M , i.e., A^ v u v — 0, n = u^N^ is the 
net baryon number density, V* = A^N U is the net flow of charge, e = u^T^Uy is 
the energy density, p + II = — | A^ u T^ u is the local isotropic pressure plus bulk 
pressure, = u v T uX A^ is the energy flow, g M = — (w/njV 11 is the heat flow, 
7r Miy = is the stress tensor, and s = u^S 11 is the entropy density. The angular 
bracket notation, representing the symmetrized spatial and traceless part of the 
tensor, is defined by = [± (A£ A U T + A^ A u a ) - \ A^ A«J A aT . The space- 

time derivative decomposes into = u^D + y^, where D = u^o^ is the comoving 
time derivative and = A^ v d v is the gradient operator. The quantity (3 = 1/T 
is the inverse temperature and a = /i/T is the chemical potential times the inverse 
temperature, w = e + p is the enthalpy density. The coefficients «,(£, n) and f3i(e,n) 
in Eq. (6) are thermodynamic integrals and given by the equation of state. 

So far is an arbitrary 4-vector and is normalized such that u^u^ = 1 and therefore 
u^dyU^ = 0. It can be taken to be the particle flow velocity. This is the Eckart frame. 
In this frame = 0. Alternatively, it can be taken to be the energy flow velocity. 
This is known as the Landau frame. In this frame = 0. The Landau or energy 
frame variables ctj and $ are related to the Eckart or particle frame variables di and 
Pi, respectively, by 



aii - a.i = /3i - /3i = w 1 , and f3 = f3 , /3 2 = #2 ■ (7) 



3 Transport Equations 



In this section we present the transport equations for the independent dissipative 
fluxes II, q 11 , 7r^. We shall consider the Landau frame where = and thus 
g M = —(w/n)V^. We shall use the Muller-Israel-Stewart second-order theory for 
non-ideal fluids [8,10], which is based on extended irreversible thermodynamics. 
Causal non-ideal fluid dynamics rests on two assumptions: (1) The dissipative flows 
(heat flow and viscous pressures) are considered as independent variables, hence, 
the entropy function depends not only on the standard variables, charge density and 
energy density, but also on these dissipative flows. (2) In equilibrium, the entropy 
function assumes a maximum. Moreover, its flow depends on all dissipative flows 
and its production rate is positive semi-definite. As a consequence the heat flow q^, 
the bulk pressure II and the traceless shear viscous tensor -n^ v obey the evolution 
equations [10] 
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r n fl + n = -C V, u» - l -CTm, 

+in,V^9" , (8) 
r q ^ r + q" = *T - a") + \kT\»3 v (^) 

- V V, ^ - V V M n + r^g, , (9) 
A^tt^ + tt"" = 2?? V (M «"> - ^Ttt^ 

+Z 7r? V <M ^ > +2r 7r 7r^<) , (10) 

where £, and r\ denote the thermal conductivity, and the bulk and shear viscous 
coefficients, respectively. Also 



rn = (Po , r q = nTp 1 , TV = 2t7/3 2 , (11) 
Jnq = C a o , = K,Ta , l q7T = nTai , l nq = 2rjai , (12) 

are the relaxation times for the bulk pressure (r n ), the heat flux (r 9 ), and the shear 
tensor (tv) ; and the relaxation lengths for the coupling between heat flux and bulk 
pressure (lu q , l q u) and between heat flux and shear tensor (l qn , l nq ). Since the bulk 
viscosity vanishes for an ultrarelativistic ideal gas [13] and since the heat flux will 
turn out to be negligible small, we shall only need the expression for 02 . a M = w M is 
the 4- acceleration, and = A M ° A"^ d[pu a ] is the vorticity tensor. 

The standard theory of Landau is contained in these causal transport equations. To 
see this, set all relaxation time and length scales in Eqs. (8-10) to zero: 



n = -ce, (13) 
M (v^_ a ,A _^Il v , (t) , (14) 



T J w VT, 

7T^ = 277(7^ , (15) 

where 6 = V/*^ is the expansion scalar, and = yWt^ = y^u^ — | A* 1 " \/\u x 
is the shear tensor. 



4 Longitudinal boost invariance with transverse expansion 



We now consider the three-dimensional expansion with longitudinal boost invariance 
and cylindrical symmetry in transverse direction. We take the z-axis to be the beam 
direction and use cylindrical coordinates. 
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We are interested in the transverse motion. We will take advantage of the one- 
dimensional scaling solution which, under boosts along the longitudinal direction, 
is invariant. Thus we can go to a point (along the z-axis) where the equations are 
simple, for instance z — 0, and derive the radial solution. Afterwards, one can boost 
this solution along the longitudinal direction. We use the coordinate transformations, 

t = rcosh?7 , z = rsinh?7 , (16) 

where 

r = VW^, v = l^l^ z • (17) 

We take matter to be net baryon free, and thus = 0. Therefore we do not need 
to consider Eq. (1). The energy- momentum tensor in the local rest frame is 



fv> = ( £ + V±) - V ± g^ + (V z - V±) m v + (V r - V±) m»m v 

+2T rz fam^ + 2W r u^m^ + 2W z u^ , (18) 

where V r = p + r rr + IT, V z = p + t zz + IT, V± = = p + + IT and 



u"(r,77,r) = (l,0,0,0) , (19) 

P(r,77,T) = (0,0,0,l) , (20) 

m"(r,77,r) = (0,e r ,0) , (21) 

with 

e r = (cos 0, sin 0) . (22) 

The r lJ and W l are the components of and the in the local rest frame. Due 
to scaling assumption the component r rz will not contribute at z = or rj = 0. For 
the same reason W z will not contribute. 



Successive boosts with radial fluid flow velocity v r = v r (r,r) in the transverse di- 
rection and with Bjorken fluid flow velocity v z = tanhrj in the longitudinal (z-axis) 
direction give 



T v> = ( e + p ± ) u *v _ p ±g ^ + (p 2 _ <p ± ) m» + (V r - V±) m»m v 

+2r rz &m^ + 2W r u^m u) + 2W z u {v l u) , (23) 

with 



u^(r, rj, r) = 7(cosh rj, v r e r , sinhr?) , (24) 
P(r, rj, t) = (sinhr/, 0, 0, coshr?) , (25) 
m M (r, rj, r) = j(v r cosh rj, e r , v r sinh^) , (26) 
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where 7 = 1/ \Jl — v 2 . The set of three 4- vectors has the following properties: w M is 
timelike, u^u^ = 1, and the other two, / M and m M , are spacelike = m^m^ = —1. 
The three are orthogonal u^l^ = u^m^ = l^m^ = 0. The energy flux is given by 



= W r u^m v) + W z u ( n v) 



(27) 



and due to the symmetry of the problem (no energy flow in the z-direction at z — 0) 
only the radial part remains and we will simply denote it by W, i.e. W = W r . The 
shear stress tensor takes the form 

' W z cosh 2 1] — r zz r rTr ) 2 v r cos cosh i] T rr ^ 2 v r sin cosh r] W z cosh r] sinh 77 ^ 
r rr 7 2 f r cos cosh 77 n r< ^ cos 2 + IT* cos sin r rr/ y 2 v r cos sinh 77 



7T 



r rr 7 2 f r sin0cosh?7 IT* cos sin IT* sin 2 + r** r rr ^ 2 v r sin sinh 77 
^ IF r cosh 77 sinh 77 T rT/ y 2 v r cos sinh 77 r rr ^ 2 v r sin sinh 7/ IT r cosh 2 77 — T rT/ y 2 v 2 J 



where IT rz = IT r = r rr 7 2 7j 2 + r zz and IT* = IT* r = r rr 7 2 - 



(28) 



Note that in writing Eq. (28) we have omitted the r rz terms, since they will not 
contribute due to scaling at z = or at 77 = 0. For the same reason the component 
n rz will not contribute. 

The components of T^ v which will be used in the following are T 00 and T 0r . At 
77 = z = they take the form: 



T 00 = W7 2 - V r + 2Wj 2 v r , 
T 0r = W^ 2 v r + W^ 2 (v 2 + l) , 

with W = e + V r . 

The equations of motion at z = i] = reduce to 



(29) 
(30) 



<9 

_T 00 



dr 



where 



dr 



U{T^ + W)v r + V r ) 



- (T 00 +V r )S + (V r -V z )--W- , 

r r 



<9r 



-T 0r S-{V r -V 4,)--W^ , 
r r 



5 = - + - 

r r 



(31) 



(32) 



(33) 
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The velocity and local energy density are given by 



T 0r - W 



(34) 



r T 00 + V r ' 




(T 0r ) 2 - (Wf 



(35) 



For an equation of state of the form p = e/3 the energy density is determined from 



Note that in the Landau frame W = 0, and in the Eckart frame W = q. 

We now consider the transport equations for the dissipative fluxes. These arise from 
the second law of thermodynamics, Eq. (3) [8,10]. Note that, due to the symmetry 
of the problem, Eq. (9) has one independent component and Eq. (10) has two inde- 
pendent components. We shall take as equation of state p — e/3, cf. section 5, and 
this implies that we do not need to consider the equation for the bulk pressure, since 
the bulk viscosity coefficient vanishes for this equation of state [13]. Consequently, 
I! = 0, cf. Eq. (8). We also take the net charge to be zero, n — 0. It is then convenient 
to work in the Landau frame. Since n — 0, the driving force term in the heat flux 
equation (the first term on the right-hand side) vanishes, cf. Eq. (14). Although the 
remaining terms need not vanish, we checked that they are numerically small and 
thus we shall not consider the heat flux equation (9). Thus we are only left with the 
equation for the shear stress tensor, . After using the cylindrical symmetry the 
various components of Eq. (10) decouple and reduce to: 



e = — 



( 



T 00 + ^ T rr \ + 4T o 2 _ 3 ( T 0r)2 _ ^2 + gyQO^r + ^ ^2 _ ^ 




(37) 



(38) 




d 



.zz 




.zz 



T 



.zz 




(39) 



Here 
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0: 
6. 

f 



1 v r 
- + — 

t r 

d 



# + 7 



d 



d 



f , 



a ^ = 



-" + 3 e - 
r 3 

T 3 



(40) 
(41) 
(42) 

(43) 
(44) 
(45) 



The coupled system of conservation equations and transport equations is solved 
numerically using the code called LCPFCT [14]. It solves generalized continuity 
equations of the form 



I? 



1 d 



^a—l 



dr 



[r^Uv] 



~QL — 1 



{r a - 1 D 1 } + C 2 —D 2 + D 3 , (46) 



dr 



with a — 1,2,3 for Cartesian, cylindrical or spherical coordinates, respectively. The 
advantage of LCPFCT in the present application is that the additional source terms 
are included by means of the terms D\, D 2 and D%. The same set of coupled equations 
can be solved numerically using the RHLLE [15] or SHASTA [16] algorithms. 



5 Initial Conditions, Equation of state and transport coefficients 



For the initial energy density distribution in the transverse plane we use the Woods- 
Saxon parameterization 

£ ^ = e \(r-R Q) /a + l > ( 47 ) 

with Rq = 1.14A 1 / 3 fm and a = 0.54 fm. Since we consider Au+Au collisions at 
RHIC, A = 197. Here Eq is the maximum energy density for central collisions at an 
initial time To- The numerical value of Eq is determined from the equation of state 
Eo = 3aTj 4 , where a is given in Eq. (48) and Tj is the initial temperature. We will 
compare the ideal and non-ideal fluid dynamics results for similar initial conditions. 

In non-ideal fluid dynamics we also need to specify, independently, the initial condi- 
tions for the dissipative fluxes. The independent components of the shear tensor are 
initialized to r rr = = e(r)/18, t zz = — e(r)/9, where e(r) is given in Eq. (47). 
These particular values can be motivated by a microscopic model study [17]. 

For this exploratory study a simple equation of state is used, namely that of an 
ultrarelativistic gas of massless u ,d ,s quarks and gluons. Since the net charge is 
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considered to be zero, the pressure is given by p = aT 4 , and the entropy density by 
s = 4aT 3 , where 

° = ( 16+ t"')$5- < 48 > 

is a constant determined by the number of quark colors and flavors and the number 
of gluon colors. Here Nf is the number of quark flavors, taken to be 3. The only 
relaxation coefficient we need is fa which, for massless particles, is given by [10] 

& = T P ■ (49) 

The shear viscosity is given by [18,19] rj = bT 3 where 

^(^^ d+^MO - (50) 

is a constant determined by the number of quark flavors and the number of gluon 
colors. The strong coupling constant, a s , is taken to be 0.5. 

In Fig. 1 we show the profiles for the components of the shear stress tensor as a 
function of the radial coordinates at different times. 

The effect of shear stress is to reduce the longitudinal pressure (hence r zz is negative) 
and to increase the transverse pressure (hence r rr and are positive). Note that 
at any point in time the sum of the three components is zero. This reflects the 
tracelessness of the shear stress tensor. 



6 Transverse spectra and HBT radii 



To calculate single inclusive particle spectra or two-particle correlation functions in 
the fluid-dynamical framework, one assumes that a fluid element decouples from the 
fluid evolution, i.e., freezes out, as soon as the particle density drops below a certain 
critical value where the collision rate becomes too small to maintain local thermo- 
dynamic equilibrium. In our case, the particle density depends on the temperature 
only. Therefore, freeze-out happens across isotherms in the space-time diagram. 

To calculate two-particle correlation functions we use the method developed by 
Pratt [20] and others, and applied to hydrodynamics in Ref. [21]. This method is a 
generalization of the Cooper-Frye formula [22] for single-inclusive particle spectra. 
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Fig. 1. The components of the shear stress tensor in units of po = £q/3 as a function 
of r/Ro, where Rq is the initial radius of the system, at different times t = Tj + nAt, 
At = 0.025 fm/c. 



6.1 Freeze- out prescription 



Just as in ideal hydrodynamics, our fluid-dynamical model describes the time evolu- 
tion of macroscopic thermodynamic quantities which must be converted to particle 
spectra before one can compare with experimental data. To do so we use the well- 
known Cooper-Frye prescription [22] which gives the particle spectrum in terms of 
an integral over the phase-space distribution function along a so-called freeze-out 
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hypersurface S(x): 



dN dN 
E- 



d 3 p dy d 2 pj_ 



^dZ,(x)p»f(x,p), (51) 



dTj^x) is the 4- vector integral measure normal to the hypersurface and, in the case 
of ideal hydrodynamics, f(x,p) is a local equilibrium distribution function of the 
form 

r(x,p) = A [ex P (Pu^-a)-e]- 1 , (52) 

where A = g/(27r) 3 and e corresponds to the statistics of Boltzmann (e = 0), Bose 
(e = +1), or Fermi (e = —1). The degeneracy factor is g — 2J+1, where J is the spin 
of the particle. The flow velocity u p is evaluated along the freeze-out hypersurface 
E and the temperature T is calculated from the energy density on S. Longitudinal 
boost-invariance dictates freeze-out along a surface of fixed longitudinal proper time 

r /( r )- 

In general the freeze-out surface Sj is a 3-dimensional hypersurface in space-time 
which can be parametrized by three orthogonal coordinates £ M (ii, w, v). The normal 
vector on the hypersurface is determined by 

= ftT ftT-^T ' (53) 

where e puXp = —e^ v \ p is totally antisymmetric, and e puXp = +1 if (yUZ/Ap) is an even 
permutation of (0, 1, 2, 3). 



6.2 Two Particle Correlations and HBT Radii 



The two-particle correlation function for bosons can be written as [21] 
r(n ,_ 1 \kd^K,exp[ iqu ^]f(x,p)\ 2 

°2(Pl,P2j - 1 + ^ ^ , (54) 

E-\ — 5 E'} — 5 — 

'dV 2 d 3 p 2 

where K p = (pi+p^)/^ is the average 4- momentum and q p = Pi — p% is the relative 
4-momentum. The symmetry of the problem reduces the number of independent 
variables on which C2 depends. For particles emitted at midrapidity, K z = 0. Due 
to cylindrical symmetry around the z-axis we can choose the average transverse 
momentum as K± = (K, 0, 0) and the relative momenta as q out = (q ou t, 0, 0), q si d e = 
(0, q side , 0) and qi ong = (0, 0, qi ong ). 

For a fixed average transverse momentum K we define side-, out-, and long-correlation 
functions as C 2 , si de = C 2 (K, 0, q side , 0), C 2 , ou t = C 2 {K, q out , 0, 0), and C 2Aong = 
C 2 (K, 0,0, qiong), respectively. Then we define the corresponding inverse widths as 
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R side = y/hi2/q* ide , R out = yfhrt/qla and R Xong = \/\a2/q"[ ong , where q* ide is deter- 
mined from C 2 ,side(q* si de) = 3 / 2 > and similarly for q* out and q* long . 



6. 3 Dissipative corrections to the distribution function 

For a gas that departs slightly from local thermal equilibrium, we may choose (in- 
dependently at each point x) a local equilibrium distribution f eq (x,p) of the form 
(52), that is close to the actual distribution f(x,p) and set 

f(x,p) = r(x,p){l + A e mx,p)} , (55) 

where the Bose (Pauli) enhancement (blocking) factors are represented by 

A sl + e^ 1 / , (56) 

where A was defined below Eq. (52) and A eq is given by replacing / on the right 
hand side by f eq . In the calculations of single particle spectra and two particle 
correlations we shall put A eq = 1, ignoring the effects of Bose enhancement and 
Pauli blocking. This is justified because the dissipative corrections are small. 

In macroscopic dynamics we seek a truncated hydrodynamical linearized description 
of small departures from equilibrium in which only the 14 variables iV M and T^ v 
appear. The microscopic counterpart of this is a truncated description in which the 
function 

y(x,p) = kv[Av 1 f(x,p)/A] , (57) 
differs from any nearby local equilibrium value 

y eq (x,p) = a(x) - I3(x)u a (x)p^ , (58) 

by a function of momenta specifiable by 14 parameters. 

The truncated description is derived via Grad's relativistic 14-moment approxima- 
tion [9,10], or via a variational method [23], namely that y — y eq can be locally 
approximated by a quadratic function of 4-momentum, 

4>{x,p) = y-y eq ^ e{x) - e M (x)p M + e^xtyp" , (59) 

where e(x), e M (x), and e uu (x) are small. Without loss of generality e uu may be taken to 
be traceless, since its trace can be absorbed in a redefinition of e. The non-equilibrium 
distribution function (55) depends on the 14 variables a — e,(3u^ + e u , e^ v . These 
parameters are related to the dissipative fluxes. For a system which only experiences 
shear stress we just need to know [10] 

e M „ = Ctt^ , (60) 
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where the coefficient C is given by 

C = 



2w 



(61) 



7 Non-ideal versus ideal fluid dynamics 

Since we do not consider hadronization in our model studies we do not follow the 
dynamics beyond the hadronization point. We stop the dynamical evolution and 
calculate the spectra along a surface of constant energy density Ef or along given 
space-time isotherms. By truncating the dynamics before hadronization we cannot 
compare our results directly with experiment. Nevertheless, we can estimate the 
effects of entropy production in the non-ideal fluid-dynamical evolution by compar- 
ing the spectra with the ones obtained in ideal fluid dynamics. Let us note that 
an additional 7.5% [24] to 25% [25] increase in entropy is expected to arise in the 
hadronization transition, but this occurs in the non-ideal as well as in the ideal case. 
By neglecting the hadronization transition, we of course underestimate the total 
amount of produced entropy and, consequently, the final hadron multiplicity. 

We present here the fluid-dynamical expansion solutions for a Bjorken cylinder ge- 
ometry for different initial energy densities Eq and initial r . We study two scenarios. 
The first corresponds to a central Au+Au collision at RHIC, where the transverse 
radius of the hot zone is of order of 5 fm, while the time scale of local equilibration 
is roughly given by the energy loss of a parton interacting with matter, T^E/dx ~ ^-5 
fm [26]. Consequently, we choose r = 0.1i?o [21]. For the initial energy density we 
take an exemplary value of Eq = 25 GeV/fm -3 . For the second scenario we consider 
t = 1/3T which is motivated by an uncertainty principle argument [27]. In this 
case we take e = 100 GeV fm -3 , which corresponds to the maximum energy density 
reached in VNI(PCM) simulations [28]. We use the same initial conditions for the 
non-ideal and ideal fluid-dynamical calculations. In this way, the dissipative effects 
on the evolution can be observed most clearly. Since we do not consider hadroniza- 
tion we do not make an attempt to reproduce actual experimental data, which could 
require a readjustment of the initial conditions. 

In Fig. 2 we show the freeze-out time T/(r), defined by Ef or, equivalently, by Tf, as 
a function of position in the transverse plane. Due to the strong initial longitudinal 
motion, the system cools rather quickly even before the transverse rarefaction waves 
reaches the center of the cylinder. The duration of the expansion is therefore solely 
determined by the longitudinal scaling expansion. This causes the horizontal parts 
of the isotherms. The effect of longitudinal cooling is reduced for larger values of To- 
One observes that for identical initial conditions ideal fluid dynamics leads to earlier 
freeze-out. The longer freeze-out times in non-ideal fluid dynamics originate from 
the fact that the longitudinal pressure is reduced due to viscosity effects, cf. Fig. 1. 



14 



1.0 



0.5 



O.t 



T 



2.0 

(a) 

1.5 — 



i — I — i — I — r 



id. fluid 

non-id. fluid 



T f = 160 MeV 



\ T i = 500MeV|_ T. = 300 MeV 

\ \ = 0. 1 fm/c s x. = 0.5 fm/c 

\ 



J I I LI I L 



(b) 



i — I — i — I — r 



id. fluid 

non-id. fluid 



T f = 160 MeV 



i i i r i i i 



.0 0.5 1.0 1.5 0.0 0.5 1.0 1.5 2.0 



r/R, 



o 



r/R, 



o 



Fig. 2. The freeze-out time r/(r) as a function of transverse coordinate r. (a) VNI scenario, 
(b) RHIC scenario. 

In Fig. 3 we show the parton transverse momentum spectra. Due to longer freeze- 
out times in non-ideal fluid dynamics, the corresponding spectrum is enhanced at 
high p T . The reason is that non-ideal fluid dynamics creates stronger radial flow 
than ideal fluid dynamics due the increase of the transverse and the decrease of the 
longitudinal pressure, see Fig. 1. Similar results on the space-time isotherms and the 
transverse momentum spectra were found [29] by considering ideal hydrodynamics 
of a transversally thermalized system of gluons. 




1 1 

k T [GeV] k T [GeV] 



Fig. 3. Parton p± distribution for the RHIC scenario. 

In Fig. 4 we show the parton mean transverse momentum, multiplicity and transverse 
energy as a function of freeze-out temperature for the RHIC scenario. Decreasing 
freeze-out temperature corresponds to increasing freeze-out time. For both the ideal 
and non-ideal evolution, the mean transverse momentum decreases with increas- 
ing freeze-out time. The reason is that, for massless particles, the decrease of the 
temperature has a stronger influence than the increase of the radial flow. One also 
observes that the stronger radial motion in the non-ideal case increases the mean 
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transverse momentum as compared to the ideal case. Note that the initial value of 
the transverse momentum in the non-ideal case differs from that in the ideal case 
because of the dissipative corrections to the single particle distribution function. 

In ideal fluid dynamics entropy is conserved, dS/dy = constant. For massless par- 
tons, the entropy density is proportional to the number density and thus dN / dy ~ 
dS/dy. Therefore, dN/dy should be constant as a function of Tf. The 10% deviation 
from this relation observed in Fig. 4(b) for the ideal fluid case is caused partially by 
the numerical viscosity generated in the course of the evolution and partially by the 
numerical uncertainty introduced when integrating the single particle spectra over 
transverse momentum. In non-ideal fluid dynamics the entropy per unit rapidity 
dS/dy increases with time (i.e., decreases with Tf) and so does dN/dy. Longitu- 
dinal boost invariance implies that the total transverse energy per unit rapidity, 
dE±/dy, is independent of y. Due to the reduced longitudinal pressure in non-ideal 
hydrodynamics, the longitudinal work performed is also reduced and dE±/dy de- 
creases more slowly ~"+ u c — + =« — „,uu ;j„„i dy nam ics. 




0.22 0.2 
T f [GeV] 



Fig. 4. The mean transverse momentum (a), multiplicity (b) and transverse energy (c). 
The solid curves are for the ideal fluid evolution and the dashed curves are for the non-ideal 
fluid evolution. 

In Fig. 5 we show the K± dependence of HBT radii and the ratio R out /Rside- The 
consequence of the longitudinal pressure reduction is to reduce the longitudinal ra- 
dius Riong- The widths of the correlation functions correspond to the space-time 
structure of the corresponding isotherms. The effect of dissipation is to reduce R ou t 
while at the same time increasing R S ide- This leads to the reduction of the ratio 
Rout/ Rside and this works in favor of reproducing experimental data which are over- 
estimated by ideal fluid dynamics [30-33]. For discussions on the effects of first order 
viscous corrections to ideal hydrodynamical description on spectra, elliptic flow and 
HBT radii see Ref. [33]. Fig.5 is obtained with the freeze-out temperature of 110 
MeV for pions of mass 138 MeV. 
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0.2 0.4 0.6 0.8 1 
K T 

Fig. 5. HBT radii and the ratio R ou t/ Rside- The solid curves are for the ideal fluid and the 
dashed curves are for the non-ideal fluid. 

8 Conclusions and Outlook 

In this work we have compared the evolution of a non-ideal fluid of massless partons 
with that of an ideal fluid. We solved the fluid dynamical equations in a realistic 3+1 
dimensional geometry appropriate for central Au+Au collisions at RHIC energies, 
i.e., assuming cylindrical symmetry in the transverse direction and boost invariance 
in the longitudinal direction. We have used similar initial conditions in the ideal 
and non-ideal case. Since our model only contains partons and lacks a description 
of hadronization, we truncated the evolution of the system just before the onset of 
hadronization, at a (fictitious) freeze-out temperature of Tf > 160 MeV. 

We then investigated the influence of viscosity on the transverse momentum spectra, 
the parton multiplicity, the transverse energy, as well as the HBT radii. We found 
that the reduction of longitudinal pressure in the non-ideal case leads to an increase 
in radial flow. This has the consequence that, for comparable freeze-out tempera- 
tures, the mean transverse momentum is larger in the non-ideal case by ~ 100 MeV. 
Another consequence is that the longitudinal HBT radius Ri ong is reduced. Simulta- 
neously, R S ide increases while R out decreases. All this works in favour of reproducing 
experimental data, although we refrain from a detailed comparison to the latter due 
to the lack of hadronization in our model. 

Finally, dissipation produces entropy and consequently the final parton multiplicity 
increases in the non-ideal case as compared to the ideal case where it should stay 
constant on account of entropy conservation. In the ideal case longitudinal work de- 
creases the transverse energy dEr/dy with increasing freeze-out time. This decrease 
is considerably smaller in the non-ideal case due to the reduction of the longitudinal 
pressure. 
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For the future, it is important to extend our considerations to finite baryon densities, 
and to include the hadronization of partons. In order to model hadronization in 
fluid dynamics one has to consider an equation of state with a phase transition 
from hadronic matter to quark-gluon plasma. In the mixed phase the bulk viscosity 
coefficients can be large and comparable to the shear viscosity coefficient. This in 
turn will have a significant effect on the life time of the system, thereby affecting 
observables such as the HBT radii. The final hadron spectra have to be computed 
along a freeze-out surface in the hadronic phase. Only then a realistic comparison 
of the model results to experimental data is possible. 

In order to reproduce the experimental data one will be forced to adjust the ini- 
tial conditions. We expect that the latter will be very different in the non-ideal 
as compared to the ideal case. Our results for the evolution of a purely partonic 
phase indicate that a better simultaneous description of particle spectra and HBT 
radii is possible. Another important observable which should be reproduced within 
the fluid dynamical framework is the elliptic flow. This cannot be studied in the 
present cylindrically symmetric scenario; we have to generalize our study to allow 
for non-symmetric configurations in the transverse plane. 
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A Particle spectra and HBT radii from a cylindrically symmetric freeze- 
out surface 

In this Appendix we calculate the particle spectrum, and the out-, side-, and long- 
correlations in the Bjorken cylinder geometry. This is a generalization of Appendix 
B of Ref. [21] to include dissipation in particle spectra and two- particle correlation 
functions. Due to boost invariance, the surface is simply the space-time hyperbola of 
constant 17 = ^rj — zj. Note that at z — 0, 17 depends only on 77, tj = 7/(77). In 

terms of the space-time rapidity, rj = tanh~ 1 (zf /tf) we have i/(r/, 77) = T/(r/) cosh 77 
and Zf(rf, rf) = 7/(77) sinhr/. Also due to boost invariance along z, 77 cannot depend 
on 77. 

In cylindrical coordinates (r, 0, z) the freeze-out surface is parametrized as 

= (r/(r) cosh 77, r cos 0, r sin 0, r/(r) sinhr/) , (A.l) 
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with 7/(77) = tf(rf,rj = 0). The normal vector is calculated from Eq. (53) in cylin- 
drical coordinates, 

dS M = ^coshr/, — cos07p^-, — sin 0t^, — sinhr/^ r/dr 77(77) d0dr/ . (A. 2) 

The 4-momentum is parametrized in terms of the longitudinal (particle) rapidity 
y = tanh^ 1 (p z / E) and the transverse mass ttit = \p± + m 2 ] 1 / 2 as 

P 1 * = {P° , P 1 , P 2 , P 3 ) = {m ± cosh y,p ± cos 4> p ,p± sin 4> p ,m T smhy) , (A.3) 

where we have used rotational symmetry around z to choose p L = (p_i_,0, 0). The 
volume element becomes 

p M dS M = ^m_i_ cosh(7/ — y) — p± cos(0 — (f) p )^jrf(r)rdr d<p . (A.4) 

The fluid 4- velocity in the Bjorken cylinder expansion is given by Eq. (24). Conse- 
quently, 

p^Uf, = m ± 7 cosh(r? - y) - p±^v r cos(0 - <p p ) . (A.5) 

We first give here the results for an equilibrium distribution function f eq (x,p). For 
the single inclusive momentum distribution we perform the ^-integration with the 
help of Eq. (3.547.4) of Ref. [34] and the ^-integration using the formula (3.937.2) 
of Ref. [34] and then expand the Bose and Fermi distribution functions. The result 
is 



Q = E^- = AnA J2(±) n+1 1 rdr r f (r)e^ 

x|mj.Io(Q!±)Ki(zj.) - p±^h(na±)K (z ± )^ , (A.6) 

where the K u arising from the ^-integration are modified Bessel functions of the 
second kind. The upper sign is for bosons and the lower one is for fermions. The 
integration over the momentum-space angle <fi gave rise to the modified Bessel func- 
tions of the first kind l u . Here, a± = n^v r ^ and z±_ = 71^7. Note that the final 
spectrum respects the symmetries of the problem, i.e., it is azimuthally symmetric 
and boost-invariant along z. For further use, we define 

U P± ) = A^E^- . (A.7) 
dp 



Now we calculate the two-particle correlation functions in order to determine the 
HBT radii. For the numerator one has to calculate the expressions 
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J^ = A | s dS^e- x ^/ T cos[g M S^ , (A.8) 
l e 2 C! = A J^di: il K^e- KV ^ T sm[q ll W] . (A.9) 



For the side-correlation function for particles at y = we may choose K± = (K, 0, 0) 
such that pt = (E,K,q side /2,0), p% = (E, K, —q side /2, 0), K» = (E, K, 0, 0), q» = 
(0, 0, q S ide, 0), where E = [K 2 + q 2 idc /4: + m 2 ] 1//2 . Since the single inclusive spectrum 
(51) does not depend on the direction of p, the two single inclusive spectra in the 
denominator in Eq. (54) are equal. For the calculation of X^ q we perform the re- 
integration using Eq. (3.547.4) of Ref. [34]. The result is 



00 f R f 
Jf = 4vr^(±) n+1 / rdrr f 

n=l J ° 



x ^EK^zMa, b) -K-^-K (z 3 )h(a,b)j , (A.10) 



where z s = nE r y/T, a = nK^v r /T and b = q S id e rf, and the functions Iq and I n (a, b) 
(for n > 0) are defined by 



1 r n 

I = — d(fi cosh [a cos 0] cos [6 sin 0] , (A. 11) 

7T JO 

1 /-t (9™/ 

I n = — J d0cos™0cosh[acos0] cos[6sin0] = -^-^ . (A. 12) 



One can show that for the side-correlation function 2| 9 = by symmetry. Thus for 
the side-correlation function 



req\2 



KM j 

C 2 ,side = 1 + (n -eqx 2 " (A. 13) 

1^-0 ) 



For the out-correlation function the momenta are pi = (Ei, K + q out /2, 0, 0), p 2 = 
(E 2 ,K- w/2,0,0), if" = (if , if, 0,0), = (E 1 - E 2 , q ouU 0, 0), where £ 1>2 = 
[(If ± q ut/^) 2 + 1T1 2 ] 1 / 2 , K° = (Ei + E 2 )/2. The r\- and 0-integrations separate, and 
the final result for Z[ q and X 2 q reads 
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00 r R r ( 

IT = 4tt £ ( ±) n+1 I rdrr^K 



-K 



Ki(a, P)Jo(a, b) + Ki(a, P)Jo(a, b) 

dr f\ 



K. (a, (3)Ji(a, b) + K (a, (3)Jo(a, b) 



(A.14) 



2T = 47r5:(±) 



n+1 



n=l 



rdr r, K° fc^a, (3)J (a, b) + /d(a, /3) Jo(a, 6) 



/C (a, P)Ji(a, b) + K, (a, (3)Jo{a, b) 



drt 



(A.15) 



where a = nK°'j/T, [3 = {E\ — E 2 )rf, a = nK^v r /T, b = q ou tff and (see Appendix 
B of Ref. [21]) 



Ko(a,P)= / dr)cos\(3coshri}e- acoshv , 
Jo 

roo 

K.Ja,p) = / dr/ cosh™ r/cos[/5 cosh r/]e- acosh " = (-l) n 
jo 



d n JCo(a,(3) 



/Co(a/3)= / dr] sin [(3 cosh rf\e 
Jo 



-ocoshn 



n d n lC (a,P) 



KJa0)= r d??cosh n 77 sin \B cosh v ] e - acoshri = (-1) 
Jo 

1 / >7r 

Jo(a,b) = — / d0 cos [5 cos 0] cosh [a cos 0] , 

7T JO 

J n (a, b) = — [ d0 cos" cos [6 cos 01 cosh|acos0l = — ' - , 
7i Jo oa n 



7T JO 



J (a,b) = — / d0sin[5cos0] sinh[acos0] , 



1 r w 

J n (a,b) = — / d0cos n 0sin[5cos0] cosh[acos0] 

71 Jo ' 



d n Jp(a, 6) 
<9a n 



Then the out-correlation function reads 



C2,out — 1 + 



Wijx e > 2 , ± ) • 

Finally, for the long-correlation function 



(A.16) 
(A.17) 
(A.18) 

(A.19) 
(A.20) 
(A.21) 
(A.22) 

(A.23) 
(A.24) 



Jf = 47r £(±)n+l / rdrr/ 

71=1 ^ 



x |£Gi(a, c)Io(6) - K-^-G (a, c)h(b)j , 

and = by symmetry. Here = [fT 2 + m 2 ] 1 / 2 , a = nEj/T, b 
c = qiongTf and 



(A.25) 

nK'jVr/T, 
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G (a,c)= / d?7e- acosh?? cos[-6sinh?7] , 
Jo 

G n (a,c)= / dr7cosh n 77e~ acosh,? cos[-6sinh?7] = 
Jo 

Then the long-correlation function reads 

C -14- ^ 



<9a ra 



(A.26) 
(A.27) 



(A.28) 



We now discuss and give the results for the dissipative corrections to particle spectra 
and HBT radii. For shear viscosity corrections we need to compute ir fiu p fl p u . From 
the shear stress tensor (28) and the 4-momentum (A. 3) we have 



PuPv^^ = m 2 L n zr cosh 2 (?7 — y) — 2m±p±r rr cosh(r] — y) cos(0 — P ) 

+p 2 _ L U r<t ' cos 2 (0 - <j> p ) + p\t h - m 2 ± T zz . (A.29) 

Performing the r\ and <p integrations over f eq (x,p)(fi(x,p) gives 



°° r Rf 

I$ = 4TrA J2(±) n+1 / rdrry 

n=l J0 



x<mln 



2 nrz 



m ± -[K 3 (z ± ) + 3K 1 (z ± )]I (a ± ) 



p±^\K 2 (z ± ) + K (z ± )]I 1 ( a± ) 



—2m_ L p±'y 2 T rr v r 



+p 2 ± U^ 



m ± - [K 2 (z ± ) +K (z ± )]l 1 (a ± ) 
- P±^ L ^ 1 (z ± ) [l 2 (a ± ) + l (a ± )\ 

OTf 2 



m ± -K 1 (z±) [h(a±) + Io("_l)] 

- P±^\k (z ± ) [l 3 (a ± ) + 3Ix(a ± ))] 



m_LKi(2;_L)Io(a_L) 



- P±^ J -K (z±)l 1 (a±) 
art 



, (A.30) 



for the shear dissipative corrections to the particle spectra. 

For the side-correlation function the shear dissipative correction to Xf ls is given by 
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If * = 4tt £ (±) n+1 / rdrr 7 

n=l J0 



x^E 2 W z 
-2EK-f 2 T rr v. 



E^(K 3 (z s ) + 3K 1 (z a ))/ (a, 6) - A^(K 2 (z s ) + K )/i(a, 6)^ . 



^(K 2 (z s ) + K (z s ))7 1 (a,6)-XK 1 (z s )7 2 (a,6)|^ 



+ J FsT 2 n r ^ 



SKi(z 8 )/ 2 (a, 6) - XK (z s )7 3 (a, 6) 



dr f 



dr 



fi 



+ (ifV^ - E 2 t zz ) 



EK 1 (z s )I (a, b) - KK (z s )I 1 (a, b) 



drf 



, (A.31) 



while for the out-correlation function the shear dissipative corrections to Jf s and 
Zf s are given by 



Zf s = 4vr ^ (±) n+1 rdrr f 

n=l ^° 

K° (/C 3 (a, /3) Jo(«, 6) + £ 3 («, /?) Jb(a, &)) 



x K° 2 n 



9r f 



-tf (/C 2 (a, /3) Ji(a, 6) + /C 2 (a, /3) Ji(a, b)) ^ 



-2K°K^ 2 T rr v, 



7sT° (/C 2 (a, (3)Ji(a, b) + £ 2 (a, /3) Ji(a, 6 



9r f 



-X (3)J 2 (a, b) + ^(a, b)) ^ 



-K (jC (a, P)J 3 (a, b) + K (a, (3)J 3 (a, b 



d_Tf_ 

dr f 



+ (k 



2 ^<P _ X° 2 T ZZ 



K° (/d(a, (3)J (a, b) + /3) Jo(a, 6)) 

/ 



-X (/C (a, /3) Ji(a, 6) + /C (a, (3) Ji(a, 6)) ^ 



(A.32) 



and 



23 



It = 4vr]T (±r +1 / rdrr f 

n=l J0 

K° (£ 3 (a, P)J (a, b) - K 3 (a, (5)J Q {a, b 



x lK° 2 W 



-K (jC 2 (a, P)J x (a, b) - K 2 (a, 0)J x {a, b 



drj_ 
drt 



-2K°K-f 2 T rr u 



K° {£ 2 {a, /3)Ji(a, b) - JC 2 (a, /3) Ji(a, 6)) 



(/C 1 (a,/3) l 7 2 (a,6)-/C 1 (a,/?)i 2 (a,6 



drj_ 

drt 



+K 2 IT 



K° (t^a, P)J 2 (a,b) -ld(a,P)J 2 (a,b)) 
-K (/C (a, f3)J 3 (a, b) - /C (a, /3) ^(a, 6)) ^ 
+ (kV^ - i^° 2 r 2 j #°(/Ci(a,/3)JoM) - IC^a, P)J (a,b 



-K (fco{a, /3) Ji(a, 6) - /C (a, /?) Ji(a, 6 



dry 
9rf 



(A.33) 



respectively. For the long-correlation function the dissipative correction is 



Xf s = 4vr j>2 (±)" +1 f Rl rdrr f 

n=l J° 



x(E z U 



2tttz 



EG 3 (a,c)l (b) - KG 2 (a,c)h(b)-^- 



-2EK-i 2 T rr v r 



f 

EG 2 (a,cMb) - Kdia^ihib) + Iq(6))^ 

2 dr f 



+K 2 Il r * 
+ (k 2 t h - E 2 t zz ) 



EG^c^iUb) +I (6)) - ^G (a,c)i(I 3 (6) +3I 1 (6))|^ 



^(0,0)10(6) -KGofacMb)-^- 

dr f. 



. (A.34) 



The full correlation functions in side-, out-, and long-direction are given by Eqs. 
(A.13),(A.24) and (A.28), with Jf replaced by l- q +lf s . 
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